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Abstract 

The relationship between a neuron's complex inputs and its spik- 
ing output defines the neuron's coding strategy. This is frequently and 
eff'ectively modeled phenomenologically by one or more linear filters 
that extract the components of the stimulus that are relevant for trig- 
gering spikes, and a nonlinear function that relates stimulus to firing 
probability. In many sensory systems, these two components of the 
coding strategy are found to adapt to changes in the statistics of the 
inputs, in such a way as to improve information transmission. Here, 
we show for two simple neuron models how feature selectivity as cap- 
tured by the spike-triggered average depends both on the parameters 
of the model and on the statistical characteristics of the input. 
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Neuronal dynamics are characterized by nonlinearities that lead to large, 
approximately stereotyped voltage excursions, or spikes, that are the basis 
for interneuronal signaling. Capturing the relationship between inputs and 
the resulting pattern of spike outputs from a given neuron in the form of 
a reduced functional model is a focus of sensory neuroscience. In the sense 
that such a model provides a general mapping from input to output, it can 
be thought of as the neuron's "c oding strategy" . 



R everse c orrela t ion methods (IBrvant and Segundd . Il976l : Ide Boer and Kuypeii . 



19681 : ISakail . Il992l : iHunter and Korenbergi . Il986l ) provide a means to sam- 



ple the statistical characteristics of stimuli that tend to trigger spikes; in 
the simplest case, the mean, or spike-triggered average stimulus (STA), 
is t he optimal l i near k ernel for predicting the firing rate from the stimu- 



lus f lRieke et al.l . Il996l ). Using reverse correlation, one may obtain an ap- 



terms of the input features that drive the system ( 


Meister and Berrv II. 


1999; 


de Ruvter van Steveninck and Bialek. 


1988; 


Simoncelli et al.. 


2004). These 



methods may be applied not only to determine how neural systems are 
driven by external stimuli, but to extract a model for how specific patterns 
of synaptic current inputs drive single neurons. This allows one to deter- 
mine the role that a single neuron with a characteristic complement of ion 
channels plays in a circuit: the i ntegration of inputs over a certain timescale 
(ISlee et all . l2005l : ISvirskis et all . l2003l : IPrescott et aP . l2006l \. the detection of 
sudden change or hi ghly synchronous events flAbeled . Il982l : ISlee et all . l2005l : 



Svirskis et al.l. 120041). or the selection of certain frequency components in the 



20061 ). 



input (llzhikevichl . 120011 : IPrescott et al. 

Here, we will derive explicit expressions for the outcome of such a sta- 
tistical analysis applied to two simple neuron models. We have two goals. 
The first is to develop a general framework for understanding how the de- 
tails of neuronal dynamics establish or influence the features in the input 
that trigger spikes. Second, neuronal systems show adaptation to statistics, 
in the sense that the neuron's coding strategy often changes when driven 
by stimuli with different statistical properties. In the case of single neurons, 
such effects can modulate or gate the effective computation of the neuron 
according to the statistical properties of the signal or the background in puts 
f lHasenstaub et all . l2007l : lOestexhe and Pard . Il999l : IPellous et all . l2003h . To 
identify the rules governing this process, one would like to know to what ex- 
tent the observed changes may result from time-independent neuronal non- 
linearities and to what extent they must be due to changes in underlying 
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neuronal parameters. To study this, we will compute how the experimentally 
obtained features of two fixed models depend on the statistical properties of 
the stimulus, focusing on the variance of a white noise input. 
The key points of this paper are: 

• The relevant linear filter corresponding to a nonlinear spiking neuron 
model is determined by the nonlinear dynamics linearized in a manner 
consistent with the typical operating regime of the system, which is 
determined both by its dynamics and by the stimulus conditions. To 
characterize this regime, we compute the voltage probability distribu- 
tions from the Fokker-Planck interpretation of the models. 

• We then use a novel application of the technique of stochastic lineariza- 
tion to map the nonlinear models onto a set of linear models. By study- 
ing both the mapping, determined by an optimization function relating 
the linear and nonlinear models, and the related STA predictions for 
the equivalent linear system, we can delineate the roles of different 
nonlinearities on spike encoding. 

• The form of the STA is influenced both by the subthreshold (non) linear 
dynamics and the spike afterhyperpolarization. 

• Models with similar phase space topology can have STAs whose form is 
controlled by different mechanisms. A rapid-onset exponential integrate- 
and-fire model (EIF) has no significant subthreshold nonlinearity, and 
so its STA is almost completely determined by the probability cur- 
rent due to spiking. In contrast, the quadratic integrate-and-fire model 
(QIF), while superficially similar to the EIF, has an STA whos form 
is dominated by the sampling of the subthreshold nonlinearity, with 
spiking effects playing a secondary role. 



1 Models and numerical methods 



Change in the effective feature selectivity with driving varianc e has been stud- 



ied for the case of the leaky integrate-and -fire (LIF) model (jPaninski et al. 



20031 : iPaninskil . l2006bl : lYu and Led . l2003l ). In the LIF model, the dynamics 



are linear until the voltage reaches an imposed threshold after which the 
voltage is immediately reset below threshold. Thus, the LIF contains no 
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intrinsic excitability, and further, does not allow for the possibility that the 
system can cross threshold multiple times before spiking due to noisy inputs. 
This discontinous behavior with respect to spike initiation is not found in 
biological neurons. 

Two simple models with more realistic spike initiation are the qu adratic 



and exponential integrate-and - fire ra odels (QIF and EIF, respectively) (lErmentrout and Kopell 



19861 : iFourcaud-Trocme et al.l . l2003l ). Both models are similar in spirit to the 



LIE insofar as they replace the afterhyperpolarization mechanism with a dis- 
continous jump, or after-spike reset, but the point of reset in these models 
occurs at the peak of the spike instead of at the threshold voltage. This 
mitigates the effects of t he pathological b ehavior in response to noise that 
the discontinuity creates (jPaninskil . l2006bl ) by moving it away from the in- 
teresting region of spike initiation. The models are described by an equation 
of the form: 



TmV = -V + f{v) + (K - Vs) 5{V - Vs) + 



(1) 



where v denotes the membrane voltage, is the membrane time constant, 
Vs is the voltage that defines the spike height and Vr is the post-spike reset 
voltage. The input current, is a zero-mean gaussian white noise (GWN) 
process with correlation function (s(t)s(O)) = a'^Tm5{t). The delta-function 
is shorthand for the act of resetting the voltage to Vy. after it reaches Vg. All 
of the spike-generating and nonlinear subthreshold dynamics are encoded in 
f{v). Eor the two models studied here, we have: 



av 



gexp 



v-V, 



for the QIE, 
for the EIE. 



(2) 



We study parameters such that the resting potential is zero and the unstable 
fixed point is at a^^ for both models. This requires us to choose K = 
(1 + ga\n{ga)) and g <^ K- Somewhat paradoxically, despite the higher 
order nonlinearity, the choice of small g causes the exponential nonlinearity to 
turn on over a much tighter range in voltage than the quadratic nonlinearity 
of the QIE. We will see that this leads to more linear behavior of the EIE 
model below threshold. Thus the two models behave noticeably differently 
below threshold while still having the same after-spike dynamics. 
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1.1 Reverse correlation analysis 

Reverse correlation is used to determine characteristics of the stimulus that 
are correlated with neuronal response. From a long, random stimulus presen- 
tation s(t) and the resulting spike response times ti, one collects the set of 
current traces that led to a spike, s(r— tj), over an interval of time r = [0, — T] 
prior to the spike where T is chosen appropriately to capture all the stimulus 
history that is relevant to triggering the spike. The spike-triggered average 
or STA, s(r), is found by averaging these samples over i: 

1 ^ 

^(^) = ]7E^(^-^0- (3) 

i=l 



1.2 Defining spike times 

The results of reverse correlation analysis can depend on how the spike time 
is defined. Here, we will look at the STAs with a temporal resolution that 
is short compared to the average spike width. Different choices of the volt- 
age threshold used to define spikes will accordingly lead to STAs that differ 
from each other due to temporal jittering of the ensemble of spike-triggered 
trajectories. Because we want to understand how the spiking of the model 
determines the feature selected from the stimulus ensemble, we are inter- 
ested in choosing a threshold that yields an STA that best captures the role 
of the stimulus on the approach to the spike but is not sensitive to stimulus- 
driven variations in the spike itself. For both models considered here, this is 
achieved by selecting the unstable fixed point that separates the subthresh- 
old region from the spiking region in the absence of noise. Since the location 
of the unstable fixed point is a function of the mean input current and the 
quadratic form of the nonlinearity , we will call it the dynamical thr eshold 
in accordance with previous work f Izhikevich . 20001 : Hong et al. . 2007 ). For 



the zero- mean inputs considered here, the dynamical threshold is Vth = oi~^- 
Thus, we define spike times as the time of the last upward crossing of the 
dynamical threshold preceding an after-spike reset. 
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1.3 Model simulation 



In discrete-time with time step h, the nonhnear models in equation [T] were 
reahzed as: 

Vn+l = Vn + — i-Vn + fiVn)) + Cr\ —^n , (4) 
'^m V '^m 

if Vn+i > K, then 

Vn+l Preset ■ (5) 

where the ^„ are drawn from a gaussian distribution with zero mean and 
unit variance. For all figures in this paper, simulations were run with a time 
step oi h = Xm/200 until 2 x 10^ spikes were accumulated. The noise was 
generated with randn in Matlab R2007b. Parameters used in simulation: 
a = 1, = 1, K = 25, K = -0.2, ^ = ^, and K = 0.77. 



2 Numerical results 

We computed the STA numerically for a range of values of the stimulus stan- 
dard deviation a for both models. Results are shown in figure [H The STA at 
all values of a has two components: an extended feature and a sharp upward 
step at the time of the spike. For the feature, two different types of behavior 
appear. For large a, the STAs are approximately decaying exponentials for 
which, as the standard deviation increases, the decay timescale decreases and 
the amplitude increases. Thus, at larger a, the QIF and EIF models perform 
approximately linear leaky integration, where the effective leakiness depends 
on the standard deviation. For very small a, the STAs are non-monotonic, 
with the peak amplitude occurring well before the spike time. 



3 Approximate STA for finite standard devi- 
ations 

To best understand how details of the models influence the STAs, we would 
like to be able to calculate the STAs analytically. In the zero standard de- 
viation limit, the STA can be analytically calc ulated for the QIF via a large 



deviat i ons principle and path integral methods (jPaninskil . l2006al : iBadel et al. 



2008al : IWilson and Steyn- Rossi . l2008l ). but that type of analysis does not ex- 
tend to flnite a. However, path integral methods can be applied for arbitrary 
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(7 to perfectly linear models with no reset. As noted previously (IHong et al. 



20071 ). the observation that the STA is an exponential implies that the sub- 



threshold dynamics of the model are effectively linear. Since the STAs of the 
nonlinear models are roughly exponential for larger standard deviations, we 
should be able to introduce a linear approximation to the nonlinear models 
that captures the qualitative behavior of the STA and helps explain in detail 
how the STA arises from the form of the nonlinearity. 

The main idea is as follows. While it is impossible to derive complete, 
time-dependent statistical distributions for these models, we can get the 
steady-state distribution from the Fokker-Planck equation. This distribution 
gives us information about how the properties of the stimulus and spiking 
dynamics determine how the system samples its subthreshold nonlinearities. 
We will then use the steady-state distribution to map the nonlinear models 
onto linear models and thus compute an approximation to the STA. Since 
the linear model follows from the steady-state distribution, we can think of 
the linear model as describing the "time-averaged dynamics" of the nonlinear 
models. 



3.1 The steady-state distribution 



A key ingredient for understanding the behavior of the spiking models and 
for determining an analytically tractable mapping of a nonlinear model to 
a linear model is the steady-state probability distribution for the voltage 
in response to an input with given statistical characteristics. This proba- 
bility djstributioiijP Fokker-Planck equa- 



tion (iPaninski et al.l. l2003l: iBrunel and Latham! . l2003l : iLindner et al.l . 12003 ; 



Fourcaud-Trocme et al.l . |2003| ). which for models of the form given in equa- 



tion [T] is: 



dpN{v,t) 
dt 



d_ 

dv 



V - f{v) 
R{t) [5{v 



PN{v,t) 



2r 



(6) 



where R{t) is the time-dependent mean firing rate that needs to be deter- 
mined self-consistently in solving the equation. This is a continuity equation 
for pn {v,t) which expresses that the evolution of the distribution is driven 
by the deterministic nonlinear driving force, diffusion, and spiking. We are 
interested in the steady-state distribution, for which = and R{t) goes 
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to the mean rate R. Using standard methods (jRiskeru . Il996l ). one can show 
that the steady state distribution is: 



Pn{v) = I ' dv'eM^'''-'^(^')), (7) 

^ Jmax(v,Vr) 

where F{v) = J f{v)dv, and the mean firing rate is the normahzation con- 
stant. 

This distribution is the product of a Boltzmann factor, controlled entirely 
by the nonlinear dynamics, F{v), preceding a spike, and a spiking flux term 
which carries the dependence on the spike parameters Vr and Vg. Since we are 
mainly interested in behavior below the unstable fixed point, or dynamical 
threshold, and the models considered here have reset voltages, Vr, near the 
resting potential, the contribution of the spiking flux term does not depend 
strongly on the form of F{v), but does depend strongly on the location of 

Vr. 



3.2 Stochastic linearization 



We turn to a set of techniques known as stochastic linearization (SL) (see 
(jSochal . I2OO5I ) for an extensive review) to model and understand the behavior 
of the STA of the nonlinear models. In the SL approach, one seeks the 
parameters of a linear model that optimally capture the properties of the 
nonlinear model in a regime of interest. In our case, we are interested in 
the linear model that best captures the approach of a nonlinear model to 
threshold for a given input standard deviation, but we are unconcerned with 
the dynamics of the spike itself. Thus, we search for an equivalent linear 
model of the form: 

TmV = -k^V + C„ + S(t), (8) 

wher e we make no attempt to model the spike or the reset (IGerstner and Kistlei 
2002 ). This linear model is simply an Ornstein-Uhlenbeck process (iRiskenl . 
19961 ). and it has the associated steady-state probability distribution: 



kn 



-exp 



kn 



V — 



k(j 



(9) 



To determine the parameters of the optimal linear model, we must se- 
lect an optimization function that maps the nonlinear model onto the linear 
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model. There are no uniqu e methods fo r choosing optimization functions 



that will yield good results (ISochal . l2005l ). and different functions will gen 



erally yield different results. We focus on two optimization functions which 
give weight to different properties of the nonlinear model. 

3.3 Minimizing the Kullback-Leibler divergence 

The Kullback-Leibl er divergence {Dkt^ measures the similarity of two proba- 



bility distributions (jCover and Thomad . l2006l ). To map the nonlinear models 



to sets of linear models, we can use the Dkl to minimize the difference be- 
tween the subthreshold part of the nonhnear steady-state distribution and 
the matched linear model's steady-state distribution. The D^l for this prob- 
lem is: 

Dkl{pn\\pl) = / dv— — In- — , 10 

Zn ZnPl[v) 

rVth 

where Z^ = / Pn{v)- 

J —QO 

To find the optimal linear model with this criterion, we minimize the Dkl 
with respect to k^- and c„. Doing so yields 



2 {E [v''] - E [vY) 

k.E [v] , (11) 



kn 



where E[...] = Z^^ I pn{v) [...]. 

J —oo 

In the cr — s> limit, ko = 1 and cq = 0, corresponding to the classical 
linearization around the fixed point of a nonlinear model. These expressions 
show that minimizing the Dkl amounts to simply estimating the mean and 
variance below threshold. This criterion is only sensitive to the probability 
distribution itself and has no knowledge of the underlying dynamics. 

3.4 Minimizing the energy below threshold 

An alternative optimization criterion is to optimize the mean square error in 
the energy, or first integral of the nonlinear models, below threshold. The 
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energy of the nonlinear model is 



V 



E^--F{v), 



and so the optimization criterion for k„ and c„ is 

,2 X 2' 



E 



-{l-K)+cv-F{v) 



(12) 



where F{v) and -E [• • •] are defined as before. This criterion amounts to 
trying to match the Boltzmann part of the distributions, taking spiking into 
account only through the bias it provides to the expectation value. This 
piece is primarily sensitive to the specifics of the dynamics below threshold 
and is less sensitive to the overall shape of the distribution than the Dkl is. 
Minimizing / yields: 



k„ = 1-2 



E [v^F{v)] E [v^] - E [vF{v)] E [v^ 



E [v^] E [v^] - E [v^Y 
E [vF{v)] E [v^] - E [v^F{v)] E [v^] 



(13) 



Again, in the cr — > limit, ko = 1 and co = 0. In this case, we see that 
the optimal parameters are directly sensitive to the form of the nonlinearity 
below threshold and that the shape of the probability distribution only enters 
through the expectation values. 



3.5 The meanings of the optimization criteria 

The two optimization criteria give different weights to different roles of the 
nonlinearity. The Dkl criterion is sensitive to the net statistical distribu- 
tion below threshold, regardless of whether it comes about due to the spike 
or the subthreshold nonlinearity, whereas the energy criterion is primarily 
sensitive to the form of the subthreshold nonlinearity. Accordingly, we can 
expect that the linear model, for a given nonlinear model and input stan- 
dard deviation, found by the different criteria will be different. Specifically, 
for nonlinear models whose optimal linear equivalents are best described by 
the energy criterion, the parameters k„ and Cg- will be closely related to the 
form of the subthreshold nonlinearity but may not be very sensitive to the 
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overall details of the voltage distribution below threshold. In contrast, for 
models with linear equivalents that are best described by the criterion, 
the parameters may have essentially no relation with the subthreshold non- 
linearity, but rather describe global statistical properties set by the mean 
and variance of the voltage distribution. 



3.6 STA of the linear model 

To find the STA for the linear model and to compare to numerical simula- 
tions, we move to discrete-time by defining t = nh where n is an integer and 
h is the time step. For clarity of notation, we identify v{t) = v{nh) = Vn- 
The linear model in equation [8] is equivalently described by the backward 
transition probability distribution: 



T„ 



1 Vn-l 



(14) 

Also of use are the steady-state probability distribution, p{v), given in equa- 
tion [HI and the forward transition distribution, p(f„|f„-i): 



hk, ^ 



P iVn\Vn-l) = yl - — j PiVn-l\Vn). (15) 

Since the model is linear, the STA follows from the spike-triggered voltage, 
V, via: 

S{t) = TrJj{t) + k„v{t) - C„, (16) 

f 

Sn = {Vn - Vn-l) ^ + k^Vn-1 - C„. (17) 

h 

The spike-triggered voltages of the linear model can be found exactly with 
the following recipe. We start at the spike time, t = (ra = 0) — the first 
time for which v > Vth and ?) > 0. 

The mean voltage at the spike time, Wq, is given by: 

/oo 
(it;ot'op(^^o|spike) (18) 
-oo 

The spike-triggered voltage distribution, p(fo| spike), follows from the threshold- 
crossing condition. The probability of finding a voltage vq at the spike time 
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is given by the probability that Vo is above Vth, niultiphed by the probabihty 
that Vq was arrived at from voltages f _i that were below threshold, summed 
over all possible subthreshold values of 

p{vo\spike) = Zq^H{vo - Vth) / dv_ip{vo\v_i)p{v^i), (19) 

J — oo 

where p(f_i) is the unconditioned distribution of voltages prior to the spike 
and is given by the steady state distribution in equation [HI p(fo|f_i) is the 
forward transition distribution, H{vo — Vth) is the Heaviside function repre- 
senting the probability for vq to be above threshold, and Zq is the normal- 
ization constant. 

The mean voltage at the time immediately preceding the spike, v^i, is 
determined by averaging over all voltages below threshold that can transition 
to voltages above threshold in the next time step, and can be found from 



= / (it>_i t>_ip(t>_i|spike), (20) 

J — oo 

POO 

p(t;_i|spike) = Z''^H{Vth - V-i) dvop{v^i\vo)p{vo\spike), (21) 

JVth 

where Z is the normalization constant for this distribution. Similarly, the 
remaining Vn for n < —2 are given by: 

/oo 
dvn-. .dv-iVnp{vn\vn+i) ■ ■ . p(t;-i | spike) . (22) 
-oo 

Equation [22] is exact for a linear model but is impractical to use. With- 
out noticeable loss of accuracy for the simulations considered in this paper, 
numerous simplifications can be made. For a discussion of approximations 
to Vq and V-i, see the appendix. 

Via the central limit theorem, the Vn for n < —2 can be simplified as: 

dVnVnP{Vn\Vn+l)- (23) 

Since these are all gaussian integrals over an infinite domain, the mean value 
is the most likely value and so the remaining averages are arrived at recur- 
sively to give: 

- hCg 

Vn+1 — 

^n= , hkT fo'^ ^^-2- (24) 
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Using the fact that — <^ 1, this recursion relation can be solved in terms of 
the exponential function and gives: 



Co- / _ Cc 



ka{n + l)h 



for n < —1. 



(25) 



The STA immediately follows from equation [111 and is: 



2k„ ( v^i - exp 



krr {n+l)h 



n < -1, 
n = 0. 



(26) 



For the linear model, the STA is simply given by the exponential filter up to 
a singular piece at the spike time that arises from requiring that threshold 
be crossed from below. 



3.7 Comparison to numerics 

For the QIF, the energy criterion qualitatively captures the time constant 
at all a and correctly matches the amplitude of the STA at high a (see 
figures [2JA. and (Sp). Conversely, the Dkl criterion is much less accurate. 
While it too predicts qualitatively correct time constants, the amplitude of 
the predicted STAs is much too large (not shown). This is because the 
Dkl criterion strongly weights the effects of the after-spike reset and total 
probability mass, and thus biases the resting potential of the linear model too 
far below threshold. Thus, the STA for the QIF is primarily determined by 
the form of the subthreshold nonlinearity and is less sensitive to the escape 
from the subthreshold domain due to spiking. Figure [3]A shows how the 
optimal linear model relates to the subthreshold nonlinearity in this case. 

For the EIF, however, the energy criterion gives /Co- ~ 1 for all a. While 
this is not surprising given the effectively linear subthreshold dynamics, it 
does not agree with the numerics. The Dkl criterion, on the other hand, 
applied to the EIF leads to qualitative agreement between the linear models 
and numerics (see figures[23 and[3P). This confirms the idea that the changes 
in the STA in the EIF can only be due to the reduction in the time spent 
below threshold because of spiking, and that the small amount of nonlinearity 
below threshold for the parameters used is not relevant except at small a (see 
figure [3]B for further discussion) . Thes e changes in th e STA are analogous to 
those studied previously by Paninski ( Paninski et al.l . boosi ). 
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The singular upward step at the spike time arises from the condition that 
thresh old must be crossed from below — that v must be positive — to elicit a 
spike ( lAguera y Areas et all l2003l : iHong et all 120071 ). This "delta-function" 
component, shown in figures [Tp and[Tp, appears here so prominently because 
we have chosen the spike-defining threshold to be at a voltage for which 
the stimulus is still relevant to spiking. This step does not vanish in the 
continuous-time limit. The value of the step can be calculated approximately 
(see appendix) with good agreement with simulation data (see figure Hj). This 
mode occasionally appears in experimental STAs when the spike waveform 
is slow (R. Mease, personal communication). It is usually not seen because 
the threshold is generally drawn well into the intrinsically excitable domain 
of the voltage and so a condition on i) does not significantly constrain the 
stimulus in that situation. 



4 Discussion 



Due to nonlinearity, LN characterizations of neural systems show dependence 
on stimulus statistics, even w i thout changes in t he underlying dynamical 
parameters (ITheunissen et al.l. I2OOOI: IYu and Led. l2003l: iBorst et all 12005 ; 



Gaudrv and Reinagell . 120071 : iGill et al.l . l2008l : IWestwick and Kearnevl . 120031 ). 
In particular, by changing only the input standard deviation, the effective 
computation changes its functional form and timescale. We have explored 
the consequences of this for two reduced naturally-spiking neuron models, 
the quadratic and exponential integrate-and-fire models. In determining 
the linear filter or filters character izing the model, our work differs signif- 
icant l y from previous approa ches te erstner and Kistler . 2002 : Hong et al. 



2OO7I : Iwilson and Stevn-Rossl . [2008. : Badel et al.l . l2008a[ ) in that the point of 
linearization is not taken, as is classically done, to be the equilibrium point; 
rather, we allow the subthreshold voltage distribution to determine the op- 
timal point of linearization. This distribution carries information about the 
form of the nonlinearities, the mean firing rate, and the stimulus itself. These 
properties account for changes in the effective linear model with stimulus 
variance. We find that despite these models' superficial similarities, different 
mechanisms are primarily responsible for this form of adaptation. The key 
difference between the models is that the QIF is nonlinear below threshold, 
qualitatively corresponding to a neuron with hyperpolarizing currents that 
are activated below threshold, while the EIF is mostly linear below thresh- 
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old. Both model s have been successfully fit to neuronal data from a vari ety 
of neuron types dlzhikevichl . l2004l : iRauch et all l2003l : iBadel et all l2008bl l 

Thus, both the neuron's intrinsic properties and the statistics of the back- 
ground or of the driving stimulus ensemble determine the effective filtering 
properties of the system. This shows one means by which modulating the 
statistics of the input can effectively gate the trans mission of different types 



of input or stimulus features through the sys tem (iHasenstaub et al.l . 12007 ; 



Destexhe and Pard . 1 19991 : iFellous et al.l . l2003l ). While this analysis focused 



on very simple model neurons, the methods we describe generalize to more 
complex, higher dimensional neuronal models, although analytical solutions 
are unlikely. These simple examples give a clear insight into intrinsic modu- 
lation of feature selectivity. 



Our previous treatments of th i s problem (Aguera v Areas et al.l . 12003 ; 



Aguera y Areas and Fairhalll . l2003l : iHong et al 



20071 ) concentrated on the 



case where spikes are well-separated, so that the effects of spike history are 
explicitly separated from the role of the stimulus in determining the proba- 
bility of generating a spike. Another approach to this problem is to include 
an explicit spike-history term in the generative m odel (iGerstner and Kistlerl . 



2OO2I : iPaninski et al.l . l2004l : iPowers et al.l . |2005| ). Here, the spike history is 



incorporated into the computation of features due to the effects of the mean 
firing rate on the steady-state distribution of threshold escape and reset. 
These results underscore the difficulty in inferring information about under- 
lying biophysical parameters from the output of reverse correlation, indepen- 
dent of a consideration of the stimulus properties. 
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Appendix: Approximating the singular piece 
of the STA 

In numerical investigation, we find that there is a simple approximation for 
the average voltage at the spike time, Vq, for the range of a considered in 
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this paper and the use of the dynamical threshold for spike triggering. In our 
hands, this relationship does not seem to follow from an obvious perturbative 

calculation. For the QIF and EIF, we find to first order in a^/—: 



vo^V,,-^ + fa\l^ (27) 

where / = 0.85 is the result of a fit to the exact integral in equation [18] for 
different a, h, and r^. 

The exact integral for the average voltage at the time immediately pre- 
ceding the spike, -U-i, can also be approximated with an error of a few parts 
in a thousand. The distribution, p(t>_i| spike), can be approximated as: 

p(t;_i|spike) ^ — jy^^ (28) 

dv^ip{v^i\vQ) 



oo 



h 

uer in u \i 

show that: 



To first order in (T\ where Vq is given by equation [271 it is possible to 




V-i^V,,-'^-aJ-\\ - + f{--l\ (29) 



Taken together, this shows that the singularity in the STA at the spike 
time, which is given by goes as (y\f^- See figure [H for numerical results. 
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Figure 1: The STAs for (A) the QIF and (B) the EIF (using the L2-norm 
for easy visual comparison), triggered on Vth = a~^, for various a with the 
upward step at t = removed. Note that at small a, the STA is non- 
monotonic while, at large a, it is approximately a decaying exponential. As 
representative examples, the STA in real units is shown (C) for the QIF for 
a = 0.3 and (D) for the EIF for cr = 3 with the last time-step included. 



Figure 2: Comparison between the numerical STA (solid) and the STA pre- 
dicted by equation [26] (dashed). (A) For the QIF, numerical results (solid) 
are compared to the prediction from stochastic linearization via the Energy 
criterion (dashed). In this case, the D^l criterion predicts STA amplitudes 
that are too large (not shown) (B) For the EIF, numerical results (solid) are 
compared to the prediction from SL via the D^l criterion. In this case, the 
Energy criterion predicts /c^r ~ 1 for all a (not shown). 
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Figure 3: The left column refers to the QIF and its optimization via the 
energy criterion, and the right, the EIF via the Dxl- (A) This figure shows 
how the linear model corresponding to the energy criterion for a = 2 for the 
QIF corresponds to the full nonlinear model. In gray, we see the steady- 
state voltage distribution. The quadratic nonlinearity is shown dashed and 
the optimal linear model as determined by the energy criterion is the solid 
line. We see that the linear model in this case is closely related to the 
average slope of the quadratic nonlinearity below threshold. (B) In contrast, 
for the EIF, using the Dkl criterion, the optimal linear model does not 
correspond closely with the exponential nonlinearity. This is indicative of 
the fact that the adaptation of the STA in the EIF is due primarily to the 
spiking reset, as evinced by the STA results (see figure [2]). (C,D) The steady 
state distributions of the (QIF, EIF) models (solid lines) are compared to 
their linear model approximations (dashed). 



Figure 4: The value of the STA in the singular component at the time of the 
spike, s(0), for all cases studied in this paper. As explained in the appendix, 
the value is approximately linear in a and model- independent. 
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